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ABSTRACT 

Effects of a turbulent gas flow on the stability of a 
glow discharge are investigated. While it is known that 
turbulence affects the stability of glow discharges, the 
mechanism is not clear. Primarily, the problem lies in the 
very long characteristic times of turbulence phenomena when 
compared to the glow discharge instability times. 

A mathematical model is developed which is solved numer- 
ically for the ambipolar diffusion and an overall energy 
equation in an unsteady, cylindrical coordinate system. 

Strong perturbations of the electric field are introduced 
which disturb an otherwise stable configuration and the 
effects of turbulence on the time evolution of the perturba- 
tion are observed. 

It is shown that the modification of the ambipolar 
diffusion coefficient and the thermal conductivity is a 
reasonably sufficient model to introduce the turbulence 
effects. It is found that the charged particle density is 
effectively the sole source of heating in the gaseous environ- 
ment of the discharge. It is then shown that turbulence 
acts to suppress the temperature instability introduced by 
the discharge streamer. 
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INTRODUCTION 



I . 

If a gas such as nitrogen or argon, or a mixture of 
gases such as air, is subjected to an electrical potential, 
a sequence of phenomena take place as a result of a continu- 
ously increasing applied voltage and/or current. Initially, 
for a small voltage, a very weak current on the order of 0.1 
picoamperes develops. This current will increase slightly 
with increasing voltage until all the naturally occurring 
charge carriers are being utilized. A current saturation 
condition now prevails until the applied voltage reaches 
about 6 kv in air at which time a further voltage rise 
rapidly increases the current flow to about one picoampere. 
This is referred to as the Townsend Discharge Region and it 
is due to new charge carriers being created by collisions 
between the highly energetic electrons (fast moving) and the 
neutral (slow moving) atoms. This process can have several 
effects on the individual gas molecules. Essentially, it 
can excite them, dissociate them, and/or ionize them. This 
ionization effect is higher for lower gas pressures since 
the farther the electron travels between collisions the 
higher will be its kinetic energy which is imparted by the 
acceleration produced by the applied electric field. When 
the electron's energy is greater than that required for 
ionization of the neutral gas particles the collisions will 
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produce additional electrons which in turn gain energy, 
produce collisions, and therefore, increase the total number 
of electrons present. This process is called an electron 
avalanche . 

From this point, slight increases in voltage produce 
current increases of several orders of magnitude at which 
time the current becomes self sustaining. Reducing the 
external circuit resistance now increases the current flow 
and a steady glow is produced in the gas. This glow- 
discharge is the familiar phenomenon which occurs in neon 
lights and has other practical applications such as a 
pumping scheme for gas lasers. 

Finally, when the current is increased still further, a 
highly ionized channel is formed which collapses the glow. 

If there is sufficient current carrying capacity in the 
external circuit a steady current is established which 
represents a plasma state called an arc. 

A full understanding of the arcing mechanism is diffi- 
cult due to the extremely short formation times involved. 
However, since delaying the point of transition to arcing is 
a highly practical problem, a method for modeling it is 
needed. A blend of theory and experiment is usually re- 
quired. Such a model could then be used to study the 
effects of the arc on and by its environment. 

It has been observed that this glow-to-arc transition 
can be delayed by utilizing gas dynamic stabilization 
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techniques in flowing gas systems [Ref. 1]. The purpose of 
this work is to establish the extent of fluid dynamic 
stabilization that results from the introduction of varying 
levels of fluid dynamic turbulence into the flow. It will 
be shown that this method has beneficial effects on the 
stability of the glow discharge process and on the temperature 
distribution in the gas produced by the presence of discharge 
streamers. For further discussion and background information 
relating to ionization and the breakdown processes refer to 
any standard reference in the field or to Wallace, R.J. 

[Ref. 2]. 

An introductory skill level in atomic physics plus 
fundamental proficiencies in mechanics, electrical princi- 
ples and numerical analysis should provide the reader with 
adequate preparation for understanding this material. 

The principles discussed here have many practical engi- 
neering applications. For maximum power output of an elec- 
trically pumped gas laser, such as a carbon dioxide laser, 
the glow discharge must be held essentially at the limit of 
arcing. Since turbulent flows tend to delay this glow-to-arc 
transition, the laser can be allowed to operate at either 
the highest possible power level or with a greater margin of 
safety at lower power levels. Another application where a 
stabilized high frequency electrical spark discharge could 
be used as a means for injecting heat into a propulsive duct 
(like a pulsejet) has been suggested by C.E. Tharratt [Ref. 3]. 
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This type of a device is theoretically capable of unlimited 
thrust. It would also have a low specific fuel consumption 
and could operate outside the audible range. The possibility 
of nearly silent operation would be highly advantageous as 
an aircraft propulsion device for the future. 

Two major cases will be discussed. The first case 
concerns the behavior of a highly transient discharge in a 
non- ionized gaseous environment such as the lightning dis- 
charge. This first case concerns a pre-glow discharge type 
phenomenon since the streamer propagates in a non- ionized 
medium. However, it serves to establish the effects of 
turbulence in a simpler situation in order to form a starting 
place for this research. The second case will look at dis- 
charge streamers in a partially ionized gas which might be 
used for the pumping mechanism of a gas laser. The method 
of approach to these problems will be developed in the next 
chapter where a mathematical model will be presented. The 
governing partial differential equations for the model pre- 
sented will be solved by numerical integration techniques 
utilizing a computer. A working Fortran program is presented 
which can be used to investigate a wide variety of the aspects 
of this problem. Significant results are then presented in 
graphical format to demonstrate the stabilizing features of 
the turbulent flow. 
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II. ANALYSIS 



A. THE CONTINUITY EQUATION 

Highly ionized channels which develop at the final 
stages of a glow discharge just prior to arc breakdown are 
called streamers. The transition to arcing is then often 
due to streamer-breakdown. It should be noted here that the 
type of streamers to be analyzed are the non-self-sustaining 
discharge streamers that precede the arc, rather than the 
electron avalanche (onset streamers) type that occur earlier 
in the glow discharge phenomena. 

Streamers are seen to issue from either electrode and 
appear always to precede the breakdown. These streamers 
propagate virtually instantaneously for our purposes as an 
ultrafast ionization wave. They consist of a partially 
ionized plasma column for which the conductivity is dependent 
directly on the charged particle density within the streamer 
[Ref. 4]. It will be shown that the collisions made by the 
charged particles within these streamers are the main source 
of heating of the surrounding gas. The value of this 
charged particle density may be determined by the applica- 
tion of the following continuity equation to the breakdown 
process . 
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Equation (2.1) gives the time rate of change of the charged 
particle density, n, as a function of radial position, r, 
and time, t, for the decaying streamer. It is written in 
cylindrical coordinates due to the natural form of the 
streamer. A schematic drawing of the streamer geometry is 
depicted in Figure 2.1. The right hand side (RHS) of 
Equation (2.1) consists of three classes of terms. The first 
term represents the effect of diffusion due to the high 
density gradient within the plasma channel. D is the ambi- 
polar diffusion coefficient. This ambipolar diffusion results 
in the electrons and ions diffusing out of the original 
channel at essentially the same rate. Briefly, the lighter 
electrons having higher thermal velocities leave the plasma 
ahead of the slower, heavier ions. This produces a local 
positive space charge which retards the loss of electrons and 
accelerates the loss of the ions. The net effect is that the 
ambipolar diffusion coefficient is approximately twice the 
diffusion coefficient of the ions alone [Ref. 4] . 

Next consider terms 3 and 4 of the RHS of Equation (2.1) 
which represent two and three-body recombination respectively. 
Recombination is a de- ionization process which is then 
responsible here for the reduction in conductivity of the 
partially ionized gas. These processes occur when an electron 
and an ion collide at relatively low velocity and recombine 
to form a neutral atom. In order to conserve linear momen- 
tum a third body must be present. In two-body recombination 
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this third body is a photon and, therefore, this process is 
sometimes called radiative recombination. For three-body 
recombination the third body is usually an ion or neutral 
atom [Ref. 4]. Due to the large complexity of electron loss 
mechanisms accurate values for the recombination coefficients 
are difficult to obtain. The values used in this thesis 
will be discussed in the next subsection. 

The second term on the RHS is the convection term. This 
term is neglected since the characteristic time of convection 
(length/flow velocity) is much less than the characteristic 
recombination times associated with terms 3 and 4 [Ref. 2]. 

While the convection term in Equation (2.1) can be 
neglected, the diffusion term cannot since the diffusion 
coefficient D a can increase by about three orders of magni- 
tude between laminar and turbulent flow [Ref. 2]. This will 
cause the characteristic diffusion time to be on the same 
order of magnitude as the recombination time and some rela- 
tionship will exist between diffusion and recombination which 
will affect the charged particle density distribution, n, 
in space and time. 

Once Equation (2.1) is solved for n(r,t) the resulting 
distribution profile can be integrated to determine a 
conductance of the decaying streamer as a function of the 
streamer radius and time. This can be accomplished under 
varying turbulence level conditions by changing the effec- 
tive value of the diffusion coefficient. 
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Equation (2.1) describes the lightning-discharge type 
phenomenon in a non-ionized medium. A more interesting 
situation is the case of the streamer as it appears in an 
ionized-sustainer discharge for the pumping of a gas laser 
[Ref. 5]. For this case, Equation (2.1) may be written as. 



3n _ ^a 3__ , 3n , 

3 1 r 3r 3r ; 



2 

v . n - a^n 

l 2 



3 

a^n 



( 2 . 2 ) 



where now a new first term is added to the RHS which repre- 
sents a production of electrons and ions. The ionization 
coefficient, v^, is a function of the sustainer E/N and can 
be determined from Figure 2.2. The electric field referred 
to here is an existing axial electric field, that is, the 
"pumping field" that maintains the glow. Equations (2.1) 
and (2.2) are non-linear partial differential equations 
which will be solved numerically by the use of a finite 
difference approximation scheme. 

B. THE ENERGY EQUATION 

In addition to the continuity equation, a suitable 
energy equation must be considered in order to describe the 
temperature distribution introduced into the gas by the glow 
discharge streamer. For thermal stability to exist, a 
balance must be created such that the time rate of change of 
the internal energy is zero because the loss of heat through 
thermal conductivity and the addition of energy from Ohmic 
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Figure 2.2. Relevant Properties of Nitrogen 
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heating must balance each other. The following equation 
represents this relationship. 
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The first term on the RHS of Equation (2.3) represents 
Joule or Ohmic heating where, E, is the electric field 
strength and, , is the current density. The second term 
is the heat flux loss due to thermal conductivity and 
radiation. The total derivative. 
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3 0 

reduces to, -r—, since the model under consideration here 
0 t 

allows for conduction only in the "z" direction. It is 
further assumed that the change of the internal energy along 
the length of the streamer is small, the radiation emitted 
primarily perpendicular to the streamer centerline has a 
negligible effect and, due to the short charac teristic times 
involved, the heating occurs at an essentially constant 
volume. With these assumptions in mind. Equation (2.3) 
becomes , 



pC v H - K,2t = ° E l 



(2.5) 



where, p, is the gas density, a, is the electrical conduc- 
tivity of the plasma (recalling that j = aE^) and K is the 
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thermal conductivity. Equation (2.5) gives the temperature, 

T, as a function of the radius and time. 

The conductivity, a, is given by, 

a = e n p (2.6) 

e e 

where, e, is the elementary electronic charge in coulombs, 

n , is the number of electrons and, u , is the electron 
e e 

mobility. The mobility is also a function of E/N and is 

defined as the ratio of the electron drift velocity, v . , to 

d 

the applied external electric field. The value of is 
given in Figure 2.2 as a function of E/N. 

2 

It will be shown later that for the oE term from 

z 

Equation (2.5) , E is essentially constant and since a is 
directly proportional to n the solution of Equation (2.2) 
represents the single source of heating that drives the 
temperature distribution obtained from Equation (2.5). The 
consequence of this is that a mathematical model describing 
the temporal and spatial characteristics of the streamer can 
now be developed. Furthermore, the effects of turbulence as 
a stabilizing medium can then be studied since the thermal 
conductivity, K, and the ambipolar diffusion coefficient, 

D , are the only parameters affected by turbulence in Equations 
(2.1) , (2.2) and (2.5) . 

C. SOLVING THE EQUATIONS 

An effective means of generalizing a problem such as 
this is to non-dimens ionalize the equations. To begin then, 
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after some algebraic manipulation. Equation (2.1) can be 
written in the form, 
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where n = n/n and r = r/r . "n " and "r " are the initial 
o o o o 

values for the charged particle density and the streamer 
radius. The following additional definitions will prove 
useful in this model, 



3 



a 



2 

a.n r 
2 o o 



D 

a 



a-,n 
3 o 




( 2 . 8 ) 



(2.9) 



/v 



t E 




( 2 . 10 ) 



The parameter, 6/ is an inverse measure of the turbu- 
lence level since the ambipolar diffusion coefficient is 
smaller in laminar flow (large 8) and larger in a turbulent 
flow (small 6). If one neglects the diffusion term in 
Equation (2.1), then for two-body recombination, direct 
integration yields, 
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( 2 . 11 ) 



This ion loss formula [Ref. 2] then shows that the 
characteristic recombination time, t r , is simply l/c^r^ . 
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Incorporating the above definitions allows Equation (2.7) 
to be written in the non-dimensionalized form, 



3n 1 3 3n 

— ~ g — ( r — 

at p 3r 3r 




( 2 . 12 ) 



which lends itself to a numerical solution by substituting 
the appropriate finite difference approximation formula 
[Ref. 6] for the differentials. After some manipulation. 
Equation (2.12) can be approximated by, 






(2.13) 



where , 
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In this scheme, ' i' represents position perpendicular to 
the streamer centerline and 1 j' represents relative position 
in time after some initial condition. Equations (2.13) 
through (2.16) form the basic algorithm used by the Fortran 
program discussed in Appendix A. 



24 



The accuracy of the charged particle density obtained by 
this type of differencing scheme depends on the step sizes 

A /\ 

At and Ar . The steps must be smaller than the smallest 
characteristic time of the phenomena described. Therefore, 

A 

At is taken as 0.1 or in real time x R /10 [Ref. 6] . Primarily 
for programming convenience, Ar is also taken as 0.1. Addi- 
tionally, the value of 3 is limited by the following stability 
criterion [Ref. 2]. 



^ > ± 77 (2.17) 

Ar 

Returning to Equation (2.2), it is seen that the addition 
of the ionization term to the RHS requires a modification 
to Equation (2.16) which now becomes, 
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Equation (2.5) can be non-dimensionalized in a similar 
manner by writing it in the form, 
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where x a is defined as, 
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A check of Equation (2.20) reveals that t has the 
units of seconds and, therefore, is a characteristic time 
for the energy equation. If one now defines the following 
additional relations , 
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then, it can be shown that. Equation (2.19) becomes the 
following non-dimensional relation, 
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dt A 3 A r 8r 8r ° E Z ' 
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(2.24) 



where, T , is the initial temperature, which for this work, 
is taken as 300° Kelvin. It should be clear here that, 8 , 

r\ 

is now a non-dimensional parameter whose value depends only 
on the thermal conductivity, K. This parameter then serves 
the same function for the energy equation as does 8 for the 
continuity equation, that is, it allows the influence of 
turbulence to be introduced into the model. The value of 
the electric field strength, E^, from the RHS of Equation 
(2.24) can be estimated by assuming that the initial value 
of the externally applied field, E q , is given by [Ref. 8], 
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A reasonable value for ft/£ might be on the order of 
10,000 ft/m, for which Equation (2.25) becomes, 
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where I/E for the streamer can be shown to be, 
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and , 
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"N" is the total number of particles present in the gas at a 

pressure, P, and a temperature, T. The Boltzman constant is 

k = 1.38E-23 J/°K. The integral in Equation (2.27) is a 

measure of the conductance of the streamer and is evaluated 

in Chapter III by integrating the charged particle density 

distribution obtained from Equation (2.1). The initial 

2 

value of this integral is simply Trn o r Q . Using estimates 

of n and r to be 3.I25E18 m~ 3 and 4.0E-4 m respectively 
o o 

[Ref. 7], Equation (2.26) can be solved for E^ becoming, 
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(2.29) 
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o 



E 



z 



1 + 0.0006 



o 



o 



where the value of y^N in Equation (2.27) can be obtained 



Results to be discussed in Chapter III show that the value 
of the integral drops rapidly to a small and nearly steady 
value (on the order of 0.01) such that the denominator of 
Equation (2.29) approaches unity in very short times. 

The results of the previous paragraph then show that the 



RHS of Equation (2.24) is simply n since o/a is (en y /en y ) 

o eeoe 

2 2 

and E /E is approximately one. So now finally, Equation 



Equation (2.30) then shows that the temperature distri- 
bution due to the presence of the streamer is affected only 
by the charged particle density. 

As for the continuity equations, Equation (2.30) can be 
solved by introducing the appropriate finite difference 
formula as follows (after some rearranging of terms) , 



from reference tables to be 5.257E24 (v-m-s) 1 for nitrogen. 



o 

(2.24) can be written as, 
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(2.31) 



where , 
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(2.33) 



The Fortran program discussed in Appendix A solves 
Equation (2.31) simultaneously with Equation (2.1) or (2.2) 
depending on the case being considered by using the new 

A 

value of n calculated from each step as an input to the 
calculation of the new temperature value. The time step for 
the temperature calculation is adjusted such that the real 
time step for the two equations is equal. 



D. MODELING THE STREAMER 

The results of Chapter III show that appropriate selec- 
tion of the controlling parameters for Equations (2.2) and 
(2.30) will produce a stable system in space and time. Then 
the question is: what effect does the momentary appearance 

of a streamer have on this stable situation? If a local 
perturbation of the electric field occurs, then the ioniza- 
tion coefficient will momentarily be raised which will in 
turn raise the charged particle density in the region of the 
perturbation. Figure 2.2 provides the data for estimating 
the magnitude of this perturbation by noting the relative 
increase in E/N for a corresponding increase in v^/N. The 
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ratio of these values can then be used to determine a 
simple multiplier (labelled Z and Z ) for the n terms 
appearing on the RHS of Equations (2.2) and (2.30) to simu- 
late the effect of the streamer. Therefore, one can write, 



3 _n _ 1 9 9_tk 

3 1 8 r 3r 3r 



A A A A 

Zn - n - an 
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(2.34) 
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Z A n 



(2.35) 



where Z and Z^ are the multipliers discussed. For the 
finite difference approximations, these two additions are 
easily handled for Equations (2.2) and (2.30) respectively, 
as follows. 




an 



3 

i , j 



Z 



. n . 

l l 



(2.36) 



T. -.1 = At[-==-(A+B) + Z n . . ] + T . . 

i/3 + l 6 a A 1,3 1,3 



(2.37) 



The Fortran program incorporates this concept by utilizing 
a separate subroutine called "STREMR" which is identical to 
the subroutine "DENSTY" except for the addition of Equations 
(2.36) and (2.37). In this scheme then, the streamer is 
modeled by calling up this subroutine when desired and then 
observing the reaction of the system over time. 
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III. RESULTS AND DISCUSSION 



A. NON-IONIZED ENVIRONMENT CASE 

In this section the effects of turbulence on the 
streamer in a non-ionized environment are considered. 

Equation (2.12) applies to this situation for which Figure 
3.1 shows a typical plot of the charged particle density 
versus radial position perpendicular to the centerline in a 

fully laminar flow. Note that the abscissa and the ordinate 

/\ 

are expressed in the non-dimensional coordinates of r = r/r Q 

and n = n / n Q ' respectively. There are five curves plotted 

which show the charged particle density distribution profile 

from the initial conditions at time zero in steps as specified 

in the legend to the profile at 3000 recombination times (TR) . 

"TR" is the computer's version of "t " and will be used here 

in the text. Now, for this plot, and except for Figures 3.11 

and 3.12, for all the figures in this section, n Q = 3.125E18 
-3 3 

m and = 2.0E-12 m /sec [Ref. 7] . These values then 
equate to a recombination time of 0.16 microsecond, and so, 
1000TR is equivalent to 0.16 millisecond in real time. 

In the laminar case depicted in Figure 3.1 it is clearly 
seen that the centerline charged particle density drops off 
in time in accordance with Equation (2.11) which in non- 
dimensional form would be, 



n 



1 



1 

A 

1 + t 



(3.1) 
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For example, at 1000TR, = 1.0E-3. This decay is due 

primarily to recombination only since, in this case, 
beta = 50,000, and therefore, the diffusion term becomes very 
small. For longer times, diffusion comes into play and some 
spreading of the charge profile is seen. 

Figure 3.2 shows the effect of introducing increasing 
amounts of turbulence into the flow. This plot compares the 
different charged particle profiles as they would appear at 
0.32 msec for beta values of 50 (fully turbulent), 500, 5000 
and 50,000 (fully laminar). Note that the 2000TR curve from 
Figure 3.1 is identical to the beta = 50,000 curve in Figure 

3.2. The effect of the turbulence is clearly seen to spread 

/\ 

out the profile and to lower n on the centerline as well 
[Ref . 2] . 

As was mentioned previously, the conductance of the 
decaying streamer can be determined by integrating the 
charged particle density over the radius. This can be 
derived by "inverting" Ohm's law and introducing the conduc- 
tivity expression for a partially ionized plasma for which 
the current can now be expressed as [Ref. 4] , 



I = GV 



(f.v 



e^n A 

tVT )v 

e e 



(3.2) 



where G is the conductance, or the reciprocal of the resis- 
tance, in the familiar relation, V = IR. 'o' is the 
conductivity, 'A' the cross-sectional area and is the 
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length of the streamer (or the distance between the elec- 
trodes in Figure 2.1 for this case) . Additionally, m^ and 
v e are the electron mass and the electron collision frequency, 
respectively . 

For the decaying streamer, the radius, and therefore the 
area, is changing but the length is constant. Since all the 
other factors within the brackets of Equation (3.2) are also 
constant, 'G' can then be expressed as. 



G 



Jim v 

e e 




n 2 it r dr 



(3.3) 



where gamma is now defined as the value of the integral in 
the above equation. Gamma has the units of reciprocal 
meters, and is determined in the program, "CHARGE" (Appendix 
A), by the subroutine, "CNDUCT." Gamma is then directly 
proportional to the conductance. 

Figure 3.3 shows gamma plotted over time out to 100TR. 
Note that gamma drops in value by 90% in only about 10TR 
(1.6 microsec). After this time the value levels off to 
about 0.04 at 100TR and then beyond the range of Figure 3.3 
the value drops very slowly to about 0.001 at 10,000TR (1.6 
msec) . 

Figure 3.3 is for fully laminar flow. Introducing 
turbulence into the flow causes a very slight increase in 
gamma over all time for the complete streamer. Figure 3.4 



35 



GRMMR (1 /METER) 



REPRODUCED AT GOVERNMENT EXPEN 



GRMMR VS TIME 
IN LRMINRR FLOW 



a 




FIGURE 3.3 



36 



GflMMfl U'/METER) ><10 

22.0 23.0 24.0 25.0 26.0 



REPRODUCED AT GOVERNMENT EXPENSE 



GflMMfl VS BETA 




FIGURE 3 . 4 



37 



shows this effect at 1000TR. Figure 3.5 shows the value of 
gamma within the original streamer tube channel only. Here 
it is seen that for small beta (high turbulence) the conduc- 
tance is lower than in the laminar flow condition. 

Referring back to Figure 3.2, and with these conductance 
characteristics in mind, it seems rather apparent that one 
effect of the turbulence is to delay the possibility of a 
transition to arcing since the conductance is lower along 
the original channel. Another consideration is what happens 
to a follow-on streamer appearing after the charged particle 
density profile of the first streamer has died away for say, 
2000TR as in Figure 3.2. In a laminar flow there would be 
a distinct path of least resistance along the original 
channel. However, in the turbulent flow, the profile is 
spread out over a much larger radius and then, this "follow- 
on streamer" essentially has no favored "route" for arcing. 

Equation (2.30) showed that the charged particle density 

was the sole source of heating in this model. Figure 3.6 

shows the result of solving Equation (2.30) in a fully 

2 

laminar flow for an E/N of 1.635E-17 V-cm , at one atmosphere, 
and 300°K for an electric field strength of 4 . 0E4 V/m. Note 
that the temperature rise amounts to less than about 0.05% 
in 100 OTR and then rises very slightly over time eventually 
peaking out at about 4000TR and then falling back slowly. 

This would seem to be a logical behavior for this case since 
the charged particle density that is heating the gas is also 
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dying away. The question is then; what happens if one 
increases the electric field? 

Figures 3.7 through 3.10 show the effects of increasing 
E q from 1.0E6 V/m (a typical value of interest for lasers) 
to 8.0E6 V/m. Note that 8, changes for each of these cases. 
This is due to the fact that the electron drift velocity is 
a function of E/N, which is plotted in Figure 2.2. Since the 
drift velocity changes, the electron mobility changes which, 
in turn, produces a new conductivity value for each different 
electric field strength case. Recalling Equations (2.20) 
and (2.23), one can see that each case has a different 
characteristic time and a different 6, . For increasing E , 
x A gets smaller while 6^ gets larger. 

The main difficulty with this is that the time step for 
the energy equation solution must be adjusted to ensure that 
the continuity equation is "feeding" the charged particle 
density values into the energy equation at the correct rate. 
For the continuity equation, the non-dimensionalized time 
step is 0.1 (t /10 in real time). Therefore, the true time 
for one time step (one "DT" in the program) is simply 
DT xt . Similarly, the true time for the energy equation is 
DTA x x . It follows then that the value needed for DTA is 
simply (DTxt r )/t a . 

Returning to Figures 3.7 through 3.10, it can be clearly 
seen that increasing E q increases the temperature rise, as 
would be expected. The effects of longer times is the same 
as the case for Figure 3.6 but requires much longer computing 
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times to reach the peak and to begin falling back. The 
effect of turbulence on the cases of interest, say 4 . 0E4 
V/m to 1.0E6 V/m, is to essentially eliminate the temperature 
rise . 

One remaining question for this non-ionized case is what 

effect would increasing n^ have? Figure 3.11 shows that at 

E = 1.0E6 V/m and n increased by a factor of 10 to 3.125El9/m^ 
o o - 1 

produces a more squared off profile than that depicted in 
Figure 3.7 suggesting that the temperature dies rapidly 
outside of the immediate influence of the streamer. Figure 
3.12 shows the effect of another 10-fold increase in N. 

Here the resulting temperature rise drops off at an even 
steeper rate at 2000TR than for Figures 3.7 and 3.11. 

Recalling that x R = l/c^r^, it is clear that increasing n Q 
decreases the characteristic recombination time significantly 

A 

which, in turn, lowers n such that there is a noticeable 
reduction in the temperature profile outside of this immedi- 
ate streamer influence. 



B. IONIZED ENVIRONMENT CASE 

In this section the effects of a glow discharge in an 
environment where a steady glow exists is considered. For 
this situation, Equation (3.4) is the applicable non-dimensional 
form of the governing continuity equation. 
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The temperature effects are then determined by Equation 

A 

(2.30) where n is obtained by solving Equation (3.4). The 
first question in this case might be; what effect does 
turbulence have on the temperature profile due to this 
steady glow discharge? At sub-atmospheric pressure, it is 
known that the charged particle density profile takes on a 
Bessel function distribution [Ref. 9]. Therefore, to inves- 
tigate this situation, Equation (2.30) can be solved inde- 
pendently from Equation (3.4) by calling up the subroutine 

"DENSTY" in the program and sending it a constant, for all 

/\ 

time, value for n. The Bessel function serves as the initial 

/\ 

condition profile for the value of n. This distribution, 

A 

labeled "Nl" in the program, allows n to go to zero at r = 4.8. 
Figure 3.13 shows the results of such a case for strictly 
laminar flow. The values of the controlling parameters for 
Figures 3.13 through 3.16 are: n Q = 3.125E18 m , 

3 

(*2 = 2.0E-12 m /sec and E q = 1.0E6 V/m. It is clearly evident 
that the temperature for such a laminar case would rise to 
approximately 5400 °K (recall that T/T^ = 1 is equivalent to 
300 °K) in only about 100 0TR. Since molecular nitrogen (^J 
begins to dissociate at about 4000 °K, the situation depicted 
in Figure 3.13 is undesirable. However, introducing turbu- 
lence into the flow tends to significantly reduce this 
temperature distribution. Figure 3.14 shows that a moderate 
degree of turbulence can be selected which, essentially, 
holds the temperature profile near the laminar profile without 
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going unstable as in Figure 3.13. Figures 3.15 and 3.16 
show the effects of increasing the turbulence level still 
further. Note that in Figure 3.15, the temperature profile 
is nearly "frozen out" at about 2400°K, while Figure 3.16 
is holding steady at around 1200°K. 

This asymptotic behavior of the temperature rise in 
reaching this stable profile appears to be related to the 
characteristic time associated with Equation (2.30). Observe 
that for all three figures (3.14, 3.15 and 3.16) the tempera- 
ture rise stabilizes between about 2000TR and 3000TR. This 
effect will continue to be evident in the remaining cases to 
be discussed. The parameters for the case presently being 
considered are such that the true time of 2500TR is 0.4 msec. 
This is approximately 4.5 times the characteristic time step 
in seconds for the energy equation, which is consistent with 
the concept of a rise time in an asymptotic system. 

The final case to be discussed here is perhaps the most 

interesting. This will be the case where there is a steady 

level of ionization in the gas and a perturbation of the 

electric field is introduced to model the effects of the 

streamer. For the steady state situation at the centerline, 

B n 

=0, so one can write Equation (3.5) . 



a 3 , 3n. 

— Tr r 3?^ 
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(3.5) 



Neglecting the three-body recombination term (a^ is on the 

6 g n 

order of 1.0E-34 m /sec) , and noting that on the centerline 
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is also zero, the ionization term is then equal to a 0 n . 

z £ 

Since Figure 2.2 gives v^/N, the charged particle density 
for the steady state condition can be obtained by. 



n 



<L 



V . 

1 




( ^i) ( JL) 

N M a 2 ; 



(3.6) 



where N is given by Equation (2.28). Figure 3.17 shows this 
steady state charged particle distribution where, for this 
model, the glow discharge is contained in a region of a 
radius equal to 10 r . At r Q , the boundary condition is set 

/N 

such that n is zero at the "wall." The left hand boundary 

condition assumes that the slope of the charged particle 

density is zero. See Appendix A for details concerning how 

the program meets these boundary conditions. Note that in 

Figure 3.17, the steady state condition is reached in less 

than 100TR, and is constant at the longer times. The first 

vertical line, at r/r =1.0, is the initial condition 

o 

/\ /s 

distribution which assumes n = 1.0 out to r = 1.0 and zero 
beyond this point. 

The steady state solution for Equation (2.30) is given 
in Equation (3.7) . 



T 



- ( B A n) 





+ T 



o 



(3.7) 



The boundary conditions are T = T at r = b and — = 0 at 
1 o ~ ~ 

/s 3r 

r = 0. Figure 3.18 shows this result for = 1.47E5 V/m. 
For this case R = 0.2219, and so at the centerline r = 0 
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and substituting in the appropriate values into Equation 

A /S 

(3.7) one obtains T = 6.548 which is confirmed at r = 0 in 
Figure 3.18. Note that this steady state condition is 
reached at approximately 2500TR which is about the same time 
period required for stabilization from the previous case 
for the Bessel function initial condition profile (Figures 
3.14, 3.15 and 3.16) and for the case showing the temperature 
profile of a lone streamer in a non-ionized environment 
(Figures 3.6 through 3.12). 

Figures 3.19 through 3.28 show the results of introducing 

the electric field perturbation into the model. Each plot 

shows the situation of the streamer appearing in a steady 

state turbulent environment and lasting for increasing 

periods of time beginning at 5 recombination times and being 

incremented by 5TR up to a final value of 50TR. In real time, 

this would be equivalent to 0.6 msec to 6.0 msec. For these 

plots, E q = 1.47E5 V/m which produces an E/N of 6.0E-17 
2 

V-cm . This case was chosen for example purposes since the 
steady state temperature profile is convenient at a centerline 
value of T/T q = 6.5 and the rise due to the perturbation 
reaches a maximum of about 13.75 (4l25°K) at 50TR. Each figure 
shows five curves, 3 are solid and 2 are dashed. The solid 

A 

lines represent the initial temperature distribution (T = 1.0) , 
the steady state profile, and then, the profile of the streamer 
at the end of its stay time. The dashed lines are the decay- 
ing, or transient, profiles produced after the appearance of 
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the streamer. Note that the centerline temperature drops to 
about half of the rise value in only 20TR and has nearly 
reached steady state in another 100TR. 

Figures 3.29 and 3.30 show a case for a higher E/N value 
2 

of 7.0E-17 V-cm . Here the steady state temperature profile 
has a centerline value of about 10.0 and then a streamer "on 
time" of only 8TR sends the profile above 4000°K. Cases for 
smaller E/N values produce lower steady state temperature 
profiles. The streamer effects are similar except that 
much longer "on times" are needed to produce much of a 
temperature rise. This appears consistent in the model 
since at the lower temperature it would take a much larger 
perturbation to upset the steady state condition, while at 
the higher energy levels a comparatively smaller perturbation 
is all that is required. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



This thesis presents a mathematical model for describing 
the effects of varying levels of turbulence on glow dis- 
charge phenomena. The development is concerned only with 
the aftereffects of the streamer, rather than its propaga- 
tion. The program that is presented in Appendix A is 
intended to allow the investigation of a large variety of 
cases of interest relating to the glow discharge. The cases 
discussed in the results chapter are those that appeared to 
develop a logical sequence of thought in attempting to 
answer the question of how turbulent gas flows can and do 
affect the stability of a glow discharge and the problem of 
delaying the glow-to-arc transition. 

The results clearly show that the modelled turbulent flow 
does indeed have strong beneficial effects on the stability 
of the system. Additionally, the mathematical model success- 
fully predicts within a reasonable level of uncertainty a 
region of stability that is consistent with the normal oper- 
ating ranges of, say, an electrically pumped gas laser 
[Ref. 10]. The model suggests that the region for which the 
system maintains thermal stability is critically dependent 
on E/N. The research showed that a steady state temperature 

increase from 1200 °K to 3000 °K results from an increase in 

2 2 

E/N of only 5.0E-17 V-cm to 7.0E-17 V-cm . 
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One important question still remains however. That is; 
is there a particular level of turbulence which would have 
the most beneficial effect on the streamer-initiated break- 
down problem? This work suggests that the best situation is 
simply that the most desirable level is the highest level 
that can be maintained. However’, the non-linear nature of 
the equations suggests that some optimal level may exist 
[Ref. 2]. Experimental work is probably the best means of 
answering this question. 

The results presented here by no means exhaust the aspects 
of this problem that could be investigated. The intent of 
the program is that it can be used as a tool to analyze any 
desired case. Some logical next steps might be to expand on 
the streamer modeling aspects to investigate what happens when 
a follow-on streamer appears before the initial one has fully 
decayed back to the steady state condition. The "Z" and "Z 
arrays could be set up to have two or more streamers appear 
at any desired location. Also, the values of "Z" and "Z A " 
used -in this model are somewhat arbitrary and probably have 
to be specified for each individual E/N case according to 
experimental observation. More research into the actual 
magnitude of the electric field perturbation may allow the 
model to more accurately predict the areas of instability 
with streamer on times that are more in agreement with 
experiment [Ref. 5]. 

Mathematically, the analysis of this problem consists of 
solving the continuity equation for three cases. The first 
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is for the RHS of Equation (3.4) equal to less than zero. 

This is the non-ionized environment case. The second is for 
the RHS equal to zero. This is the steady state case which 
is what was applied to the ionized environment cases discussed. 
The third case, which was not applied in this thesis, is for 
the RHS of the continuity equation to be greater than zero. 
Investigation into this third case could possibly clarify 
some aspects of the problem. And finally, the energy equation 
may have eigenvalue solutions which might produce results of 
interest such as additional regions of stability. 
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APPENDIX A 



THE CHARGE PROGRAM 

The Fortran program "CHARGE" is intended to be a working 
tool for analyzing the glow discharge problem. It evolved 
from a previously written program which solved the non-linear 
continuity equation (2.1) [Ref. 2] . It is now arranged in a 
somewhat different form. The program consists of a main 
program section which calls up any needed combination of the 
seven subroutines which perform the desired tasks. As cur- 
rently written, it can compare four separate cases plus it 
will display the initial conditions of either the charged 
particle distribution or the temperature profile. Individual 
test cases are set up by commenting out (placing a "C" in 
column 1) those commands not desired. After the main program, 
there are two sections. The first provides commands to print 
out the data and the second provides commands to plot the 
results by utilizing the "DISSPLA" system developed by the 
Integrated Software Systems Corporation, San Diego, California. 

The program as it appears following this discussion is 
set up to produce the plot displayed in Figure 3.28. The 
main program begins by defining all arrays. For this exam- 
ple the arrays are set to 101 in the "i" direction which 
specifies the size of the array with respect to radial posi- 
tion. Since one DR, one radial step, is equal to r Q /10, 
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position "101" specifies a distance of 10 r Q perpendicular 
from the centerline of the streamer. The "j" size is always 
kept at 2 due to the computing scheme used in the subroutines. 
After setting the array size, set the variable named 
"RADIUS" equal to the same value as the "i" coordinate above. 
In this case, RADIUS = 101. This determines how far in the 
radial direction that the differentiation is to procede. 

Next one must assign the desired values for BETA and BETAA 
depending on the case at hand. The BMIN, BMAX, BAMIN and 
BAMAX parameters apply when it is desired to transition the 
calculation from a laminar to a turbulent environment. This 
feature is designed to attempt to ensure that the calculation 
accounts for the fact that the streamer initially appears so 
quickly that the flow looks laminar to the streamer no matter 
what the level of turbulence is. Then, after a period of 
time, the charged particle distribution left behind by the 
streamer has existed for sufficient time to now "see" the 
turbulent effects. The time required for this transition is 
estimated by noting that the characteristic diffusion time 
for laminar flow is approximately 0.1 msec [Ref. 2]. Depend- 
ing on the case at hand, one DT is on the order of a few 
microseconds. Then, 100 DT (10TR) is on the order of a few 
tenths of a millisecond. Therefore, for most cases, a 
transition of 100 DT is used. 

Returning to the example, one procedes by selecting the 
desired initial conditions. In this example, all of the 
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"Nl" values are commented out since they apply to the Bessel 
function initial condition profile. Four series of subrou- 
tine calls are provided then to allow running of four 
separate cases. Each case assigns its respective results to 
the "VALX" (X = 0,1,2,3,4...) array for plotting later. 

The initial conditions are assigned to VALO or VALTO. 

Then the first subroutine call will see the initial distribu- 
tion values of n while each successive call will see the 
results of the previous call. If it is desired to recover 
the initial conditions , call either "NZERO" or "TZERO" as 
desired and the current values will then be reset to time 
zero . 

For this example now, "VALl" is the result of two calls. 
The first is for "TRNSTY" to transition from laminar to 
turbulent conditions with time set to 101, then "DENSTY" is 
called and runs for 23,900 DT. The result is the steady 
state temperature distribution in turbulent flow at 2500TR. 
VAL2 is then the result of calling up the subroutine " STREMR" 
and running it for 500 DT with BETA and BETAA set to the 
laminar values. VAL3 then results from calling up "TRNSTY" 
for 100 DT and "DENSTY" for an additional 100 DT at the 
turbulent values. Finally, VAL4 simply runs "DENSTY" for 
another 1000 DT in turbulent conditions to check the decay- 
ing temperature profile. Of course, any combination of these 
subroutine calls and time values can be selected to view the 
situation as desired. 
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The results can then be printed out and/or plotted. To 
become familiar with the "DISSPLA" plot system, refer to 
DISSPLA Pocket Guide (current with version 9.0) by ISSCO 
Graphics, 10505 Sorrento Valley Rd , San Diego, California. 

Some pointers of potential errors should be clarified. 
When selecting the desired heading lines, ensure that bhe 
last entry in the call to "HEADIN' 1 correctly indicates the 
total number of heading lines. If no heading appears, this 
is probably the reason. Secondly, ensure that the third 
entry in the call to "CURVE" is equal to the "i" value on 
the arrays. This is the number of points (101 here) to be 
plotted. And finally, the second entry in the call to 
"LEGEND" is the number of lines of legend needed, 5, in this 
example . 

The subroutines "DENSTY," "STREMR" and "TRNSTY " all 
contain essentially the same algorithm with modifications to 
perform their particular function. The basic idea here is 
that the calculation procedes from the first row of data 
and fills the second row with updated values. Then the first 
row of old data is replaced by the second row of new data by 
"SWITCH" and the loop increments one DT and procedes. With 
this scheme the program can be run as long as desired without 
using large values of memory storage. 

The boundary conditions for the problems are met by the 
routine. The left boundary condition assumes that the slope 
is zero across the streamer centerline. So here, the new 
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value at n(2,2) is assigned to n(l,2) before the switch. 

The right boundary assumes that n is zero at the wall, that 
is station 101. This is met by limiting the maximum radius 
value to be one less than "RADIUS" (called "TEMRAD"). The 
routine procedes to "TEMRAD," then the last point seen by 
the calculation algorithm is zero. Similarly ,. the tempera- 
ture of the right boundary is set to 1.0. 

To compute the conductance parameter, gamma, that would 
result from any particular n distribution, call "CNDUCT" 
after assigning the desired n distribution to some "VALX." 
The single value "NN" will be returned, which is gamma. 

It is hoped that the program is sufficiently self- 
explanatory. All variable, constant and subroutine names 
are defined and documentation comments are provided. 
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APPENDIX B 



HP-41 PROGRAM 



In the course of developing the results of Chapter III, 
it is necessary to perform a large number of simple but 
tedious calculations. The following HP-41CV program proved 
helpful. It should be easily adaptable to any small pro- 
grammable calculator. 

2 3 

The program prompts the user for E/N (V-cm ) , v^/N (cm /s) , 
-3 3 

N (cm ), a .2 (m /s) and v^ (cm/s). These data are available 
from Figure 2.2. 

-3 3 

The program then displays E (V/m) , N (m ) , (m /s) , 

DTA, BETAL (laminar), BETAT (turbulent), BETAAL (laminar), 
BETAAT (turbulent) and DT in true time in microsecs. The 
minimum BETA values computed are consistent with the 
stability criterion. Equation (2.17). It is good practice 
to add about 1% to these values to ensure a smooth running 
calculation . 
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